N-body simulations for coupled dark energy: halo mass function and density profiles 
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We present the results of a series of iV-body simulations in cosmologies where dark matter (DM) 
is coupled to dark energy (DE), so easing the cosmic coincidence problem. The dark-dark coupling 
introduces two novel effects in N-body dynamics: (i) DM particle masses vary with time; (ii) gravity 
between DM particles is ruled by a constant G* , greater than Newton's constant G, holding in other 
2-body interactions. As a consequence, baryons and DM particle distributions develop a large scale 
bias. Here we investigate DE models with Ratra-Peebles (RP) potentials; the dark-dark coupling 
is set in a parametric range compatible with observations, for as concern background and linear 
perturbation properties. We study the halo mass function, the halo density profile and the behavior 
of the non-linear bias. We find that non-linear dynamics puts additional constraints to the coupling 
parameter. They mostly arise from density profiles, that we find to yield higher concentrations, in 
coupled RP models, with respect to (uncoupled) dynamical DE cosmologies. Such enhancement, 
although being a strong effect in some coupling parameter range, is just a minor change for smaller 
but significant values of the coupling parameter. With these further restrictions, coupled DE models 
■ with RP potential are consistent with non-linear observables. 
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The nature of Dark Energy (DE) is one of the main puzzles of cosmology. DE was first required by SNIa data 
, pj, but a flat Universe with OT ~ 0.3, h ~ 0.7 and flbh 2 ~ 0.02 is also favored by CMB and LSS observations Q 
(ft m ,b. matter, baryon density parameters; h: Hubble parameter in units of 100 km/s/Mpc; CMB: cosmic microwave 
background radiation; LSS: large scale structure). 

DE could be a false vacuum; from the expression of its stress-energy tensor, TjS — Ag^ (A is a positive constant 
and is the metric tensor), one immediately appreciates that its pressure and energy density (jpde and pde) have 
ratio w = —1. False vacuum, however, requires a severe fine tuning at the end of the EW transition. Otherwise, DE 
I could be a scalar field <j> self-interacting through a potential V (</)). Then 
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p DE = 4> 2 /2a 2 + V(<t>), p DE = ft /2a 2 -V(<j>) (1) 



' (here dots indicate differentiation in respect to the conformal time r). However, if the kinetic component Pk{4>) = 
P" 1 , (j) 2 /2a 2 approaches the potential component V(<fr), pde vanishes and the scalar field (j> behaves as CDM (Cold Dark 
Matter). This happens, e.g., in the well known case of axion DM. But, even for lower pk, if 1/2 < pt/V < 1, it is 
— 1/3 < w < and the model, at most, approaches an open CDM behavior. The relevant domain is attained when 
Pk/V <C 1/2, although keeping a state parameter w > — 1, Then, <j> approaches a ACDM behavior and is currently 
dubbed dynamical DE 0| 0| 0| 

The conceptual contiguity between DM and dynamical DE suggests that they may not be disjoint entities. If so, 
one could hopefully ease the cosmic coincidence problem, i.e. that DM and DE densities, after being different by 



orders of magnitude for most of the cosmic history, approach equal values only in today's Universe. The simplest way 

glial, 



to deal with this idea amounts to admit an interaction between DM and DE [g| . In a number of papers 
it has been shown that, at the linear level, this causes no apparent conflict with LSS or CBR data. We shall refer 
to models where DM and dynamical DE interact as coupled DE models. Other models that propose a direct link 
between DM and DE invoke a unified model (e.g. 9J) or condensation mechanisms (lcj . 

An open question, then, concerns the emergence of non-linear structures in these models and how easy it is to fit 
observed LSS data with predictable features. In this work we deal with this problem by using N-body simulations of 
coupled DE models for DE self-interacting through a Ratra-Peebles 0|[ RP hereafter] potential: 

1/(0) = A 4+a /cj) a . (2) 

Once the exponent a and the DE density parameter ^Ide are fixed, the energy scale A is set. This self-interaction 
allows w -C —1/3, if A is sufficiently low. 

Let us outline soon that our non-linear treatment sets precise constraints to coupled DE parameters. A wide 
parameter space however remains where, apparently, these models fit LSS data as well as those with (uncoupled) 



dynamical DE, although coupled DE, with RP potential, allows no improvement of such fit. The above motivations 
for coupled DE P3 however remain and, altogether, the non-linear test is successfully passed. N-body simulations 
of models with (uncoupled) dynamical DE were recently performed |iJ|- Here we follow the same pattern of 

and use the program ART ^| providing, first of all, the fair dependence of the matter density parameter fl m 
on the scale factor a. To our knowledge, this is the first time an TV-body simulation with species-dependent scalar 
gravity is carried out. Our conclusions are based on simulations of a variety of models with different RP slopes a 
and coupling parameters (3. Let us list them soon: first of all, we test two a slopes: 0.143 and 2. The latter value 
approaches the greatest value for which agreement with CMB observations is granted This is the range of RP 
models which are most distant from ACDM. We explored also a wide set of (3 couplings, ranging from 0.05 to 0.25. All 
simulations were performed starting from the same random numbers and, for the sake of comparison, we also run a 
ACDM simulation starting from such random numbers. The other parameters were set to values chosen in agreement 
with recent CMB experiments Q: tt c h 2 = 0.15, Q, b h 2 = 0.01, h = 0.7 (fi c is the (cold) DM density parameter). All 
models are normalized so that a% — 0.75 today, to match both CMB data and the observed cluster abundance [2lj | . 
Further details on the simulation performed are listed in Table 1. 



Model 


a 




Box size (h : Mpc) 


# of particle 


Mass res. (/i _1 M ) 


Force res. (h 1 kpc) 


RPi 


0.143 


0.05 


80 


128 3 


2.0xl0 10 


5 


RP 2 


0.143 


0.10 


80 


128 3 


2.0xl0 10 


5 


RP 3 


0.143 


0.15 


80 


128 3 


2.0xl0 10 


5 


RP 4 


0.143 


0.2 


80 


128 3 


2.0xl0 10 


1.2 


RP 5 


0.143 


0.25 


80 


128 3 


2.0xl0 10 


1.2 


RP 6 


2.0 


0.15 


80 


128 3 


2.0xl0 10 


5 


ACDM 








80 


128 3 


2.0xl0 10 


5 


Table I. 



For all these models we also run an high-resolution simulation of a single halo, with a mass resolution of 2.5 x 
1O 9 /i -1 M and a force resoltuion of 1.2h~ 1 kpc. 

In the next section we discuss the linear and post-linear aspects of coupled DE, explaining, in particular, how 
f2 m (a) is obtained and used. In section 3 we focus on the newtonian regime for coupled DE models and describe the 
different gravitation of baryons and DM. In Section 4, we implement these prescriptions in the numerical code, so 
explaining which further modifications ART needs, to deal with coupled DE. Section 5 is then devoted to illustrating 
the results while, in Section 6, we draw our conclusions. 



II. BACKGROUND EXPANSION IN MODELS WITH COUPLED DARK ENERGY 

Quite in general, energy density and pressure of each component, in models with dynamical DE, are obtainable from 
the stress-energy tensors TjfiJ b ' r ^\ (for CDM, baryons, radiation, and DE, respectively; radiation includes neutrinos). 
General covariance requires that the sum T M1/ of these four tensors fulfills the continuity equation 

= (3) 

and, although this is true if all tensors fulfill it separately, such requirement is not necessary; e.g., when we deal with 
fluctuations before hydrogen recombination, only the sum of baryon and e.m. radiation tensors fulfills it. In a similar 
way, if DE and DM interact, we can have that 

rW" = ^|V (C) ^ (4) 

(here is the trace of the CDM stress-energy tensor), and the sum of DM and DE fulfills eq. ©. No analogous 
interaction should involve baryons, for which we assume that T^]/ 1 — 0; in fact, experimental and observational 
constraints restrict an hypotetical DE-baryon coupling to (3b < 0.01 [ijj; also radiation cannot be involved, as the 
trace of its stress-energy tensor vanishes because of its equation of state. The particular form of the coupling <@J 
reduces to Brans-Dicke scalar gravity upon a conformal transformation (see, e.g. ^JEEiHl)- 



We assume a flat conformal background metric ds 2 — a 2 {—dT 2 + 5ijdx' l dx : >) — 1, ..,3), so that the background 
continuity equations read 

4> + 2H<p + a 2 V^ = y/l6TrG/3/3a 2 p c , (5) 

p c + mp c = -v/l67rG/3/3p c 0, (6) 

Pb + m Pb = 0, (7) 

Pr +AH P r = 0. (8) 

being H = a/a the conformal Hubble parameter. The dimensionless constant (3 2 can be seen as the ratio of the 
DE-DM interaction with respect to gravity In these models, after equivalence, the world passes through three 

different expansion regimes, denominated: (i) 0-MDE, (ii) tracking phase (@), and a final (hi) global attractor. 

Immediately after equivalence, the world enters a </>-MDE stage (^-matter dominated expansion, see |fi]), not far 
from a matter dominated expansion (MDE), although one should not neglect a small correction, proportional to (3, 
due to the kinetical term pf. (</>). V((f>), instead, is negligible and, accordingly, the correction is indipendent of the 
potential shape. By solving the Friedmann equation we find that 

a 0^/(6+4/3^ (9) 

i.e. that the scale factor grows more slowly than in a pure MDE (in Section V.B, we shall see that, on the contrary, 
the perturbation growth is enhanced, during this stage). During the (/>-MDE stage, V{cf)) gradually increases and, 
eventually, approaches and exceeds Pk(4>)] then, the world enters the tracking phase, whose details depend on the 
potential shape; for most potentials, such phase ends up into a global attractor, when DE density overwhelms DM 
and any other densities. 

Along the expansion history, the scaling of p c (DM density) is modified with respect to the uncoupled case and 
reads 

„ = ^ e -^(Mo). (10) 



here the subscript o indicates value at the present time r Q , while we take a Q = l. Meanwhile the baryon density grows 
usual. 

In the next section we shall see that these behaviors strongly affect the fluctuation growth, even in the newtonian 
regime. Fig. Q](left panel) shows the different trends of the density parameters for two different values of the coupling 
parameter (3. The three stages of the background evolution are clearly visible. For sake of completness we also report 
in Fig. (right panel) the a dependence of the state parameter w for different value of (3. Notice that, for RP models 
with a low a, the value of the state parameter has been quite close to —1, since a redshift ~ 3-4. 

All along these stages, however, according to the Friedman equations, the following relation holds (T^ |: 



dt _ x an b (a) 

da= x{a) = H ° (U) 

Therefore, once the a dependence of the baryon density parameter is given, all derivatives in respect to time can be 
easily converted into derivatives in respect to the scale factor. This relation will be used in the implementation of the 
ART program. 



III. DYNAMICS IN THE NEWTONIAN REGIME 



Let us now consider density fluctuations and discuss their evolution. First of all, the conformal metric must be 
modified to take into account the local gravitational fields and, in the absence of anisotropic stresses, reads 

ds 2 = a 2 [-(l + 2$)dr 2 + (1 - 2<P)d Xi dx 3 }, 

<I> being the gravitational potential. Let us use a = log a as independent variable, instead of r. Differentiation in 
respect to a will be indicated by '. As usual, fluctuations are expanded in Fourier components; let us then consider 
a component of wavenumber k and define A = H/k. 

Baryons and CDM will be considered as fluids; fluctuations will then be characterized by density fluctuations 



5 c ' b = Sp Ctb /p Ctb 



(12) 




Figure 1: Left panel: density parameters for radiation, baryons, DM and DE. Just after radiation equivalence, a DE plateau 
occurs, due to the dark-dark coupling. In the text this plateau is denominated 0-MDE stage. Right panel: evolution of the 
state parameter w for different value of f3 

( c ' b stand for CDM and bayons) and velocity fields v^' b = dxi/dr, from which we build the scalar variables 

e^ = i^f. (is) 

Scalar field fluctuations (S(j)) will then be described by the variable 



*=Vi^ (14) 

Dealing with fluctuation evolution well after recombination, we shall neglect their radiation component. 

Taking into account density inhomogeneities, the equations © yield the dependence of S c ' b and 9 c ' b on the scale 
factor a. The Friedman equation shall then be used to obtain the dependence on time. We shall omit here the general 
form of these equations (see fH), and will write them in the Newtonian limit, i.e. for small scales (in comparison with 
the horizon scale ~ H -1 ) and small velocities (in respect to c). 

The former condition tells us to consider the lowest order terms in A; in this limit, the gravitational potential fulfills 
equations that can be written in a simpler form by defining the function /(</>) according to 



V{tf>)=A exp [Vl67rG/3/(M] , (15) 

as well as the functions 

A =4 + /, / 2 = ^| + 2/ + /l ; (16) 

notice that this is no restriction on the potential shape. Then 

$ = -^\ 2 {n b S b + n c S c + 6Xip + 2Xip' -2Y 2 f 1 ip) (17) 
=3X(p-$, (18) 

while the scalar field fulfills the equation 

<p" + (2 + — )<f? + \- 2 <p - UXip + + 2Y 2 (f 2 <p - A$) = (5VL C {8 C + 2$). (19) 



These expressions have been simplified by using the variables X 2 = 8irGpk((j)) a 2 /3H 2 (kinetic density parameter) 
and Y 2 = 8irGV(4>) a 2 /3H 2 (potential densiity parameter) . It follows that if the DE kinetic (or potential) energy 
density gives a substantial contribution to the expansion source, X (or Y) is 0(1). 



In the Newtonian limit, however, we also neglect the derivatives of ip, averaging out the rapid oscillations of if 
and the potential term f 2 Y 2 (p; this actually requires that A << {fzY) -1 (remind that Y is 0(1))- Furthermore, in 
eq. iflfll) . the metric potential which is proportional to A 2 , can also be neglected. Accordingly, eq. ifTTIl becomes 

<P = -'^\ 2 (n b 5 b + n c s c ), (20) 

which is the usual Poisson equation, while eq. iflfljl simplifies into 

\~ 2 (p ~ pn c 5 c . (21) 

In the same way (see [8]), we obtain the continuity equations 

5 C " = -6 c '(l + ^-2px) + ^(l + ^)n c S c + ^{l b 5 b , (22) 

5 b " = -6 b '(l + ^) + ~(n c 6 c + n b 5 h ), (23) 
ri I 

and the Euler equations 

e cl = -e c (i + ^ - 2f3x) - ~(i + jp 2 )n c 6 c - ^n b s b , (24) 
e b ' = -e b (i + ^)-^(n c s c + n b s b ). (25) 

From the latter ones, we can derive the acceleration of a single non-relativistic CDM or baryon particle (mass m c ,&), 
assuming that it lays in the empty space, at distance r from the origin, where either a CDM particle of mass M c or 
a baryon particle of mass M b are set. 

In fact, owing to eq. fRjf. normalizing the scalar field so that its present value 4> = 0, and assuming that the 
density of the particle is much larger than the background density, it shall be 



n b s 



Pent 3H 2 a 
b _ PM b - Pb _ 87rG(5(0) 
Pcrit 3H 2 a 



S being the Dirac distribution. Then, reminding that the divergence V^' 6 — 9 c,b TL, and using the ordinary (not 
conformal) time variable, instead of a, eq. H24I1 yields 



v„i _ -bo. - W v,< - _ *e*m, (26) 

where dots yield differentiation in respect to ordinary time, 

G* = G(l + ^/3 2 ) (27) 

and H = d\oga/dt. 

We can integrate this equation, taking into account that the acceleration is radial, as the attracting particle lays in 
the origin. It will then be 

a ^ ■ , f , d(x 2 v) A , . 
d x V • v = 4-7T / dx^- — - = Attx v, 
J dx 

i) being the modulus of the (radial) acceleration (in the second term x = |x|). Accordingly, for a CDM particle, the 
desired expression of the radial acceleration reads 



y-(i- Wf ^.- 8,tf - f T ¥,, -«*, ( 28 ) 



(n is a unit vector in the radial direction; r = ax) which has various peculiarities and ought to be suitably commented. 
To this aim it is important to compare it with the radial acceleration 



. GM c e-V^^ GM b 



v° = -Hv° • n - - ^ (29) 



of a baryon particle. In the expression f28l . notice first the velocity term. This is a peculiar acceleration that exist 
even in the absence of particles displaying their attraction; its presence means that CDM matter is not expanding 
in a comoving way, due to the extra gravity it feels. Accordingly, its particles do not follow geodesies, because their 
mass changes in time, and their ordinary (not comoving) linear momentum obeys the equation 



G*M c e-VW I P* GM b 
Pc = 



Baryon particles, instead, safely follow geodesies, although feeling that CDM particle masses are varying. 

Let us conclude this section by summarising its specific findings: (i) The mass assigned to CDM particles does 

vary in time, being m c = m e~ V^3 while baryon particles do keep a constant mass, (ii) When interacting 
between them, CDM particles feel an effective gravitational constant G* = G(l + 4/3 2 /3); any other particle-particle 
interaction occurs with the ordinary newtonian interaction constant G. 



IV. METHODS 



Particle mass variations and different interaction constants ought to be taken into account in performing N-body 
simulations. They will be based on the Adaptive Refiniment Tree code (ART) ^| that has been suitably modified 
to deal with coupled DE models. The ART code starts with a uniform grid, which covers the whole computational 
box. This grid defines the lowest (zeroth) level of resolution of the simulation. The standard Particle-Mesh algorithms 
are used to compute density and gravitational potential on the zero-level mesh. The ART code reaches high force 
resolution by refining all high density regions using an automated refinement algorithm. The refiniments are recursive: 
the refined regions can also be refined, each subsequent refinement meshes of different resolution, size, and geometry 
covering regions of interests. Because each individual cubic cell can be refined, the shape of the refinement mesh can 
be arbitrary and match effectively the geometry of the region of interest. 

The criterion for refinement is the local density of particles: if the number of particles ina a mesh cell (as estimated 
by the Cloud-In-Cell method) exceeds the level of n t hresh, the cell is split ("refined") in 8 cell of the next refinement 
level. The refinement threshold depends on the refinement level. For the zero's level it is n t hresh = 2. For the higher 
level it is set to n t hresh — 4. The code uses the expansion parameter a as the variable tiome. During the integration, 
spatial refinement is accompanied by temporal refinement. Namely, each level of refinement, I, is integrated with 
its own time step Aai = Aa /2 l , where Aa = 3 • 10~ 3 is the global time step of the zeroth refinement level. This 
variable time stepping is very important for the accuracy of the results. As the force resolution increases, more steps 
are needed to integrate the trajectories accurately. Extensive tests of the code and comparison with other numerical 
N-body codes can be found in |2fj . 

Let us now describe the three main modifications we made to handle coupled DE. A first step amounts to distinguish 
between baryons and DM particles, which feel different gravitational forces. Therefore, the potential on the grid is 
to be calculated twice, so to fix the different forces that baryon and DM particles feel. All particles act on baryons 
through the usual gravitational constant G, which sets also the action of baryons on DM particles. DM particles 
instead, act on DM particles through a different interaction constant G* = G(l + 4/3 2 /3). The gravitational force is 
then computed through the usual FFT approach. 

A second step amounts to take into account that the effective mass of DM particles is time varying. Aside of the 
acceleration due to gravitation, each DM particle will therefore undergo an extra acceleration 2/3X. Besides of these 
two changes, peculiar of coupled DE models, we ought to take into account the right relation between a and t, as 
shown in eq. QT$ , where x( a ) = dt/da is given. By solving the background equations, in a suitable file we provide 
X(a) in ~ 200 scale factor values a*, that we than interpolate. 

The models listed in Table 1 were first simulated in a 80/i _1 Mpc box. We then selected the same halo in all 
simulations and magnified it. The low-resolution simulation, performed with a force resolution of 15 ft. -1 Mpc and 
a mass resolution ~ 2 • 10 10 ft _1 M , allowed us to evaluate the halo mass function. The high-resulution simulation, 
performed with a force resolution of ~ 1.2 ft. -1 kpc and a mass resolution of 2.54 • 10 h~ Mq, magnified a sphere with 
a radius of 5 ft, Mpc, centered on the halo, allowing us to compare halo profiles down to a radius ~ 5 /i _1 kpc. 

Besides of the above points, we could also test the non-linear evolution of the bias between the amplitudes of 
inhomogeneities in baryons and DM. Such bias is one of the most peculiar features of coupled DE models and we shall 
describe how non-linearity modifies it. 
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Figure 2: Mass function at z = for a = 2 and a = 0.143. For a = 0.143 we report three curves, for different values of /3, 
They are all practically indistinguishable and are well fitted the approximation of Jenkins et al (2001). 



V. RESULTS 
A. Mass function 

We identify halos in simulations by using a SO algorithm, that we shall now describe in more detail. As first step, 
candidate halos are located by a FoF procedure, with linking length A = U x d (d is the average particle separation) 
and keeping groups with more than Nf particles (U and Nf fixed herebelow). We then perform two operations: (i) we 
find the point, CV, where the gravitational potential, due to the group of particle, has a minimum; (ii) we determine 
the radius r of a sphere, centered in Cw where the density contrast is A„ (we use the virial density contrast found 
in the absence of dark-dark coupling [Isll^). Taking all particles within r we perform again the operations (i) and 
(ii). The procedure is iterated until we converge onto a stable particle set. The set is discarded if, at some stage, we 
have less than Nf particles. If a particle is a potential member of two groups it is assigned to the more massive one. 
In this work we use U — 0.2 and take Nf so to have a mass threshold 5.0 • 10 12 ft. _1 M Q . 

Fig. shows the mass function for isolated halos for models with different values of a and j3. Let us remind that 
the simulations have the same initial phases and the same value as — 0.75. Thus, the differences between models are 
only due to different couplings or w(t). Remarkably, at z = the mass-functions are practically indistinguishable: a 
mass-function has no "memory" of the past evolution. The mass-function obtained in this way is well fitted by the 
approximation provided by 23] for ACDM models (long dashed line in Fig.^J. 



B. Linear and non linear bias 

From eqs. 12212.3(1 . the linear evolution of the density perturbations can be easily worked out (in some cases ^| 
this can be done analytically). In Fig. we show S c - b as a function of the scale factor a. 

As a consequence of these dynamical equations, S c develops a bias with respect to S b , due to the extra gravity felt 
by DM. At the present epoch, this bias, found in the linear theory, is well fitted by the following empirical expression: 

»fofl~f ~ 1 + 0.015^ + 2.1^ ™ 



Both this expression and Fig. 0] show that the bias b depends on (3, while its dependence on a is very weak. Using 
the high resolution clusters we can test the behavior of the bias in the highly non-linear regime. To do so, we define 




Figure 3: DM and baryons linear perturbation growth for two different values of f3. The dependence on a is weak and could 
not be appreciated in this plot. 




Figure 4: Linear bias as a function of (3 for three values of a. Notice the very weak dependence on a. 



the integrated bias B as: 

Pb(< r) - p b 



B{< r) 



Pb Pc{< r) - p c 



where p c (< r) and Pb(< r) are calculated inside a radius r from the halo center and we use a hat ( " ) to denote 
average densities. In order to avoid problems with force resolution, the central zone (r < 10/i _1 kpc) of the halo is not 
used. In Fig. we show B(< r) for the same halo, in cosmologies characterized by different coupling parameters j3, 
keeping all other parameters equal. Fig.0shows that non-linearity significantly enhances the expected bias; however, 
at larges scales, we recover the theoretical linear value, as provided by eq. (|3fil) . 

The scale dependence of bias can also be appreciated from Fig.El where power spectra for baryons and DM, worked 
out from simulations at z = 0, are shown. 



C. Density profiles 

Let us remind again that all simulations are started from the same random numbers. Therefore, it comes as no 
surprise that they yield similar world pictures. In the ACDM simulation, we selected a halo, whose virial radius 
r v = 812/i _1 kpc encloses a mass M v = 6.45 • 10 13 /i _1 M Q . Similar halos, located in the same place, are set in all 




Figure 5: Behaviour of the integrated bias B for /3 = 0.15 and for (3 = 0.25. Notice that B tends to the predicted linear bias 
(dashed horizontal lines) at large scales. 




Figure 6: Power spectra for DM and baryon particles evaluated from the simulations averaging over 60 random observers. The 
increase of the bias at small scales appears clearly. 



other models considered. We then run new simulations of all models in Table 1, magnifying the region centered on 
this halo. To do so, short waves were first added to the initial perturbation spectrum in all simulations. 
In ACDM, the halo profile is accurately fitted by a NFW expression |24l25|: 

P(r) = 6* 

with a scale radius r s — 0.249 h~ 1 Mpc (here p cr is the critical density and 5* is a parameter which sets the halo 
density contrast). 

When the same halo is magnified in coupled DE models, we find model dependent behaviors towards the halo center. 
The essential restrictions to coupled DE models, arising from the non-linear treatment, derive from these behaviors. 
However, in spite of such model dependence in the central areas, the outer parts of halos (R > 100 ft. _1 Mpc) show 
discrepancies that, from 100/i _1 kpc to 700/i _1 kpc, never exceed ~ 10%. 

Let us now discuss the substantial model dependence found in the central region (R < 100 /i -1 kpc). It was already 
known that halos are denser in dynamical DE than in ACDM although the density enhancement is fairly small 
and hardly exceeds ~ 40%. Higher density means smaller r s . The coupled DE simulations we perfomed show that 
the dark-dark coupling tends to enhance such effect. In Fig. we overlap the profiles of the DM components of all 
our models, starting from ACDM (lower curve), up to a RP model with coupling parameter j3 — 0.25 (upper curve). 
The values of r s change from ~ 0.25 fc -1 Mpc (ACDM) to ~ 0.022 fcr x Mpc (/? = 0.25). The dependence of r s on (5 is 
plotted in Fig. 03 
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Figure 7: DM density profile for four different coupling values and for ACDM. 




Figure 8: Scale radius of a NFW profile as a function of the coupling parameter (3. 
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Figure 9: Upper panel: DM and baryons density profiles (respectively upper curve and lower curve) for (3 = 0.15 and for 
P — 0.25. Lower panel: once rescaled taking into account the different values fib and f2 c , there is no discrepancy between DM 
and baryons. 
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Figure 10: DM and baryons density profiles for a = 2. and a = 0.143 (here f3 = 0.15). 

In order to make sure that this effect was not related to some peculiarity of the halo selected, we magnified two 
other halos of a simulation with (3 — 0.25. Here we found even lower values for the scale radius r s (0.0105 h^ 1 Mpc 
and 0.0103 h~ 1 Mpc respectively). 

As a matter of fact, the profiles found = 0.25 or 0.2 can be fitted by a single power law: 

P( r ) 7 

oc r ' 

Per 

in the whole dynamical range, i.e., from r = 1.0/i -1 Mpc down to r = 0.005/i _1 Mpc (resolution limit), with a value of 
7 ~ -2.30. 

An analysis of Fig. |H1 shows that, only for j3 as low as ~ 0.1, r s attains half the value for ACDM. Accordingly, we 
may consider viable coupled RP models only with j3 < 0.1. 

Simulations distinguish baryons from DM particles, as already discussed in the bias subsection. This allows us to 
draw separate density profiles. They are shown in Fig. El for two different coupled DE models (upper panel). No 
apparent discrepancy between DM and baryon profiles can be seen: they overlap fairly well, once we rescale them 
taking into account the different values of and fl c (Fig El lower panel) . 

In Fig. El we compare the profiles of the same halo, with two different values of a (2.0 and 0.143), but with the 
same coupling (J3 = 0.15). The profiles overlap very well both for DM particles (upper curves) and for baryons (lower 
curves). We conclude that the slope of profiles depends very weekly on a. 

VI. CONCLUSIONS 

After finding that coupled DE models are consistent with those observables whose behavior can be predicted at the 
linear level 0.0. 0|. a test of their non-linear behavior had to be carried on. An optimistic hope was that coupled DE 
models helped to solve some of the long standing contradictions between observations and theoretical or numerical 
predictions (e.g. (13]) ■ In particular, one could hope to find halo profiles whose shape is not NFW (if this is still a 
problem) with a slope distribution closer to the observed ones for low surface brightness galaxies [2ft l2^ | and spiral 
galaxies [2^|. From this point of view, coupled DE with a RP potential leads to modest results. Very high coupling 
levels, instead of producing a flatter core, yield profiles still farther from observations. In all cases, the problem with 
concentration distribution is not solved, just as when making recourse to models with (uncoupled) dynamical DE. 

It should be also reminded that the RP potential considered here is characterized by very low a values. The scale 
factor dependence of the state parameter w, for such values of a, is shown in Fig. ^ and approaches -1 already at 
redshifts ~ 3-4. As far as the state parameter is concerned, therefore, these models are quite close to ACDM and, in 
a sense, could be considered a variant of ACDM models, for which DE is coupled to DM. 

In spite of the lack of improvement for what concerns slopes, N-body simulations lead to really significant results. 
First of all, the parameter space for coupled DE models is restricted to couplings (3 < 0.1, however leaving a wide 
room for significant couplings. Apart of the question of profiles, the halo mass function has been tested and found 
consistent with other DE models and with observations. Its evolution has been predicted and can be tested against 
future data. From this point of view, therefore, coupled DE passed the non-linear test. 




Performing N-body simulations was also important to test the evolution of the bias, between baryon and DM 
fluctuations, which is one of the main characteristics of coupled DE models at the linear level. Here we showed that 
the bias still exists, and actually increases, in the non-linear structures and studied its evolution. These results are 
open to an accurate comparison with data, that shall be deepened elsewhere. 

Let us however conclude these comments with a further observation. Coupled DE apparently leads to higher halo 
concentration essentially because of the evolution of the mass of DM particles and of the coupling constant between 
them. In the simulations we run, such mass depends on time and gradually decreases, as is predicted by coupled DE 
theories at the Newtonian approximation level. Accordingly, each DM particle mass was greater than today, in the 
past. Its gravity was therefore stronger. This is the reason why, although normalizing all models to the same as at 
z = 0, we produce more concentrated halos: the forces which bound them were stronger in the past than today. 

After appreciating this point, one can tentatively propose a way out for the halo concentration problem: a coupled 
DE model leading to DM particle masses which increase in time. This increase should also be fast enough to beat 
the higher gravity constant binding DM particels. This takes us back to the selection of a suitable self-interaction 
potential V(4>), which has no immediate obvious solution. This problem shall be therefore deepened in future work. 
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